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pH . Abstract 
• 1— I ' 

A new algorithm is presented to find exact traveling wave solutions of differential- 
difference equations in terms of tanh functions. For systems with parameters, the 
algorithm determines the conditions on the parameters so that the equations might 
admit polynomial solutions in tanh. 

\ Examples illustrate the key steps of the algorithm. Parallels are drawn through 

' discussion and example to the tanh-method for partial differential equations. 

' The new algorithm is implemented in Mathematica. The package DDESpecialSo- 

■ lutions.m can be used to automatically compute traveling wave solutions of nonlin- 

ear polynomial differential-difference equations. Use of the package, implementation 
issues, scope, and limitations of the software are addressed. 



> 



> 



Program summary 



Title of program: DDESpecialSolutions.m 
^ I Catalogue identifier (supplied by the Publisher): 

d ' Distribution format (supplied by the Program Library): 

Computers: Created using a PC, but can be run on UNIX and Apple machines 
Operating systems under which the program has been tested: Windows 2000 and XP 
Programming language used: Mathematica 
Memory required to execute with typical data: 9 MB 
Number of processors used: 1 
Has the code been vectorised or parallelized?: No 

Number of bytes in distributed program, including test data, etc. : 104 761 

Nature of physical problem: The program computes exact solutions to differential- 
difference equations in terms of the tanh function. Such solutions describe particle 
vibrations in lattices, currents in electrical networks, pulses in biological chains, etc. 
Method of solution: After the differential-difference equation is placed in a travel- 
ing frame of reference, the coefficients of a candidate polynomial solution in tanh 
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are solved for. The resulting solution is tested by substitution into the original 

differential-difference equation. 

Restrictions on the complexity of the program: The system of differential-difference 
equations must be polynomial. Solutions are polynomial in tanh . 
Typical running time: The average run time of 16 cases (such as Toda, Volterra, 
and Ablowitz-Ladik lattices) is 0.228 seconds with a standard deviation of 0.165 
seconds on a 2.4GHz Pentium 4 with 512 MB RAM running Mathematica 4.1. The 
running time may vary considerably, depending on the complexity of the problem. 

Key words: Exact solutions, traveling wave solutions, differential-difference 
equations, semi-discrete lattices, tanh-mcthod 
PACS: 02.70.Wz; 02.30.Ik; 02.30.Jr; 02.90.-Fp 



1 Introduction 



Since the work of Fermi, Pasta, and Ulam in the 1950s [14], differential- 
difference equations (DDEs) have been the focus of many nonlinear studies (for 
references see e.g. [16,39]). There is renewed interest in DDEs, which can be 
used to model such physical phenomena as particle vibrations in lattices, cur- 
rents in electrical networks, pulses in biological chains, etc. Unlike difference 
equations which are fully discretized, DDEs are semi-discretized with some (or 
all) of their spacial variables discretized while time is usually kept continuous. 
DDEs also play an important role in numerical simulations of nonlinear partial 
differential equations (PDEs), queuing problems, and discretizations in solid 
state and quantum physics. 

There is a vast body of work on DDEs, including investigations of integrability 
criteria, the computation of densities, generalized and master symmetries, and 
recursion operators [15]. Notable is the work by Levi and colleagues [20,21], 
Yamilov [41,42] and co-workers [4,9,10,29,30,38], where the classification of 
DDEs (into canonical forms), integrability tests, and connections between in- 
tegrable PDEs and DDEs are analyzed in detail. To a large extent, the clas- 
sification and integrability testing of discrete equations parallels continuous 
equations (for reviews and references consult [5,25,26,31]). 

A wealth of information about integrable DDEs can be found in papers by 
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Suris [32.33,34,35,36] and his book [37] in progress. Suris and others have 
shown that many lattices are closely related to the celebrated Toda lattice 
[40], its relativistic counterpart due to Ruijsenaars [28], the KvM lattice [18], 
and the two-component Volterra system [29,30]. 

Recently, variants of the tanh-method have been successfully applied to many 
nonlinear polynomial systems of PDEs in any number of independent variables 
[6,11,12,13,27]. Baldwin et al. [7] implemented the tanh- method (and also 
sech, cn, and sn-methods) in Mathematica. Liu and Li implemented the tanh- 
method [22,23] in Maple. 

While there has been considerable work done on finding exact solutions to 
PDEs, as far as we could verify, little work is being done to symbolically com- 
pute exact solutions of DDEs. In this paper we present an adaptation of the 
tanh-method to solve nonlinear polynomial differential-difference equations, 
which to our knowledge is novel. Our algorithm applies to semi-discrete lat- 
tices and allows one to find closed form solutions that are polynomial in tanh. 

Although the tanh method is easy to apply, it leads to fairly cumbersome 
algebra without guarantee that a tanh solution can be found [6]. We therefore 
present [8] the fully automated software package, DDESpecialSolutions.m in 
Mathematica, which implements the algorithm. Without intervention from 
the user, our software computes traveling wave solutions as polynomials in 
Tn = tanh^n, where = J2f=i di^i + Silj Cj^j + ^- The continuous variables 
Xi and the discrete variables rij are combined with constants q and dj, and 5 
is the phase constant. 

For systems of DDEs involving constant parameters (denoted by lower-case 
Greek letters), the software automatically determines the conditions on the 

parameters so that the given equations might admit polynomial solutions in- 
volving tanh. Obviously, since sech ,^ = 1 — tanh ^ our code can find solutions 
in even powers of sech . But our code can not find solutions involving odd pow- 
ers of sech. 

The paper is organized as follows: Section 2 gives a discussion of the tanh- 
method for PDEs with a worked example of a two-dimensional nonlinear PDE. 
Since the tanh-method for semi-discrete lattices so closely parallels the tanh- 
method for evolutions equations, we believe that a discussion of the continuous 
case increases the transparency of the more complicated semi-discrete case. 
In Section 3 we present our tanh-method for DDEs with a worked example 
of a two-dimensional Toda equation followed by almost a dozen additional 
examples. A description of our Mathematica package, DDESpecialSolutions.m, 
is given in Section 4. We discuss our results and draw some conclusions in 
Section 5. Appendix A gives the full input and output of a test case. 



3 



2 Traveling Wave Solutions of PDEs 



In this section we discuss the tanh-method as it apphes to systems of P poly- 
nomial differential equations, 

A(u(x),u'(x),u"(x),---,u(™)(x)) = 0, (1) 

where the dependent variable u has P components Ui, the independent variable 
X has components Xj, and u*^'")(x) denotes the collection of mixed derivative 
terms of order m. We assume that any arbitrary coefficients parameterizing 
the system are strictly positive and denoted by lower-case Greek letters. To 
simplify the notation in the examples, we will write u,v,w,... instead of 
Ui, U2, U3, . . . and x,y,t, . . . instead of Xi, X2, X3, etc. 



2.1 Algorithm of the tanh-method for PDEs 



Step 1 (Transform the PDE into a nonlinear ODE). We seek solutions in the 
traveling frame of reference, 

N 

e = E9^. + '^' (2) 

where the components Cj of the wave vector x and the phase S are constants. 

In the tanh method, we seek polynomial solutions expressible in hyperbolic 
tangent, T — tanh ^. Based on the identity cosh^ ^ — sinh^ ^ — 

tanh' ^ = sech^ ^ = 1 - tanh^ ^, (3) 
tanh" C = -2 tanh ^ + 2 tanh^ ^, etc. (4) 

Therefore, the first and consequently all higher-order derivatives are polyno- 
mial in T. Thus, repeatedly applying the chain rule, 

d» dT d» ^ d 



cA^-n-^ (5) 



dxj dxjdidT ' dT 

transforms (1) into a coupled system of nonlinear ODEs, 

r(r,u(r),u'(T),...) = o, (6) 

where Ui{T) corresponds to Mi(x). 

Step 2 (Determine the degree of the polynomial solutions). Since we seek 
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polynomial solutions 

Mi 

U,{T) = Y^a,jT\ (7) 
the leading exponents Mj must be determined before the aij can be computed. 

Substituting Ui{T) into (6), the coefficients of every power of T in every equa- 
tion must vanish. In particular, the highest degree terms must vanish. Since 
the highest degree terms only depend on T^' in (7), it suffices to substitute 
UiiT) = T^* into (6). In the resulting polynomial system P(T) = 0, equating 
every two possible highest exponents in every component Pi gives a linear 
system for determining the Mj. The linear system is then solved for Mj. 

If one or more exponents Mj remain undetermined, we assign a strictly positive 
integer value to the free Mj, so that every equation in (6) has at least two 
different terms with equal highest exponents in T. 

Step 3 (Derive the algebraic system for the coefficients a^). To generate the 

system for the unknown coefficients aij and wave parameters Cj, substitute (7) 
into (6) and set the coefficients of to zero. The resulting nonlinear algebraic 
system for the unknown aij is parameterized by the wave parameters Cj and 
the parameters in (1), if any. 

Step 4 (Solve the nonlinear parameterized algebraic system). The most diffi- 
cult step of the method is analyzing and solving the nonlinear algebraic system. 
To solve the system we designed a customized, yet powerful, nonlinear solver. 

The nonhnear algebraic system is solved with the following assumptions: 

• all parameters (lower-case Greek letters) in (1) are strictly positive. (Van- 
ishing parameters may change the exponents Mj in Step 2). 

• the coefficients of the highest power terms (ojMi, i = 1, • ' ' ^P) iii (7) are 
all nonzero (for consistency with Step 2); and, 

• all Cj are nonzero (demanded by the nature of the solutions we seek) . 

The algebraic system is solved recursively, starting with the simplest equation, 
and continually back-substituting solutions. This process is repeated until the 
system is completely solved. 

To guide the recursive process, we designed functions to (i) factor, split, and 
simplify the equations; (ii) sort the equations according to their complexity; 
(iii) solve the equations for sorted unknowns; (iv) substitute solutions into the 
remaining equations; and (v) collect the solution branches and constraints. 

This strategy is similar to what one would do by hand. If there are numerous 
parameters in the system or if it is of high degree, there is no guarantee that 
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our solver will return a suitable result, let alone a complete result. 



Step 5 (Build and test sohitions). Substituting the solutions from Step 4 
into (7) and reversing Step 1, one obtains the explicit solutions in the original 
variables. It is prudent to test the solutions by substituting them into (1). 



2.2 Application of the tanh-method to a PDE 
Consider the (2 + 1) dispersive long wave system [24] , 

f Uyt + VxX + U^Uy + UU^y = 0, 

\vt + U^ + U^xy + UxV + UVx = 0, 

which is related to the Eckhaus system. 

Applying chain rule (5) repeatedly to (8), we get the coupled ODEs, 

2c2T(c3 + CiU)U' + ciC2(T2 - l)f/'2 + 2clTV' 

+ (T2 - 1)(C2(C3 + c,U)U" + cjV") = 0, 
ci(l + 2(3T2 - l)ciC2 + V)U' + (C3 + ciU)V' 
+clc2{T^ - l){<oTU" + (T" - l)U"'] 



(8) 



(9) 



0, 



where T — tanh(cix + C2y + Cai + 5). 



To compute the degree of the polynomial solution(s), substitute U{T) = T^^^ 
and V{T) = i^to (9) and pull off the exponents of T (see Table 1). 
Removing non-dominant exponents and equating possible highest exponents, 



Term 


Exponents of T 


Uyt 


{Ml -2, Mi} 




{M2 -2,M2} 


UxUyj UUxy 


{2Mi,2Mi - 2} 



Term 


Exponents of T 


Vt 


{M2 - 1} 


Ux 


{Mi-1} 


Uxxy 


{Mi-3,Mi-l,Mi + l} 


UxV, UVx 


{Ml + M2 - 1} 



Table 1 

The exponents of T in (9) after substituting Ui{T) = T^' 

we find Ml = M2 (from Uyt and Vxx) or 2Mi = M2 (from UxUy or uUxy and 
Vxx) from the first equation. Then, from the second equation, we find M2 — 2 
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Ml 


= 1, 


f/(T) 


— aiQ + fliiT, 


M2 


= 2, 


V{T) 


= + a.2iT + 022^^, 


Ml 


= 2, 


U{T) 


= aio + aiiT + ai2T^, 


M2 


= 2, 


V{T) 


= a2o + 021^ + 022^^- 



(from UxV or uVx and Uxxy)- This gives us two branches, 

(10) 
(11) 

For the first branch, substituting (10) into (9) and equating the coefficients of 
to zero gives 

ci(2a22Ci + aiiC2) = 0, 
ci(2a22Ci + a^iC2) = 0, 
a2iCi + aioaiiCiC2 + aiiC2C3 = 0, 
aiiCi(a22 + 2C1C2) = 0, 
oiiCi + aiia2oCi + aioa2iCi - 2aiiclc2 + 02103 = 0, 
aiia2iCi + aioa22Ci + 02203 = 0. 



[12) 



Similarly, for the second branch, substituting (11) into (9) and setting the 
coefficients of to zero yields 



aiiai2CiC2 = 0, 

ai2CiC2 = 0, 

c^2iCi + O10O11C1C2 — 3aiiai2CiC2 + aiiC2C3 = 0, 
2a22Ci + a\^CiC2 + 2aioai2CiC2 - 2a\2CiC2 + 2ai2C2C3 = 0, 

2a22Ci + aiiCiC2 + 2aioai2CiC2 + 2ai2C2C3 = 0, (13) 
Ci(ai2a2i + 01x022 + 2aiiCiC2) = 0, 
ai2Ci(a22 + 6C1C2) = 0, 
oiiCi + aiia2oCi + aioa2iCi - 2auclc2 + 02103 = 0, 
«i2Ci + ai2a2oCi + aiia2iCi + aioa22Ci — 8ai2CiC2 + a22C3 = 0. 

For the first branch, solving (12) under the assumption that an, 022, Ci, C2, and 
C3 are nonzero, we find 

Oio = , On = ±ici, 

Cl (14) 

O20 — 2C\C2 — 1, 021 = 0, 022 = ~2CiC2. 

Substituting (14) into (10) and returning to u{x,y,t) and v{x,y,t), we get 

C3 

u(x, y,t) = ± 2ci tanh(cia; + C2y + c^t + 6), 

Cl (15) 

v{x, y, t) — 2ciC2 — 1 — 2ciC2 tanh^(cia; + C2y + c^t + S). 
In the second branch, the equation a\2C\C2 = in (13) is inconsistent with our 
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assumption that 042, 022, Ci, C2, and C3 are nonzero. This branch does not yield 
a solution to (8). 

Substituting (15) into (8), we verify that our solution does indeed satisfy the 
original system. 



3 Tanh method for nonlinear DDEs 

The tanh-method can be adapted to solve nonlinear polynomial DDEs. Apart 
from slight, yet important modifications, the steps mirror those in Section 2. 

Given is a system of M polynomial DDEs, 
A (u„+pi (x) , u„+p2 (x) , . . . , u„+p, (x) , 

Un+pi (x) , u'^+p^ (x) , . . . , u'„+p^ (x) , . . . , 

u!:ipi(x),u!:|p.(x),...,ua.pjx)) = 0, (16) 

where the dependent variable u has M components Ui, the continuous variable 
X has N components Xi, the discrete variable n has Q components rij, the k 
shift vectors p,j G Z*^, and u'^''^(x) denotes the collection of mixed derivative 
terms of order r. We assume that any arbitrary coefficients that parameterize 
the system are strictly positive and denoted by lower-case Greek letters. 

To simplify notation in the examples, we use dependent variables u,v,w, . . . 
instead of Ui, U2, "Us, • • • , continuous independent variables ,t, y,t, . . . instead of 
Xi, X2, x^, . . . , and lattice points n,m, . . . , instead of ni,n2, etc. For example, 
the two-component Volterra equation [34] , 

Un = Un{Vn - Vn-l) 
Vn = Vn(Un+l - Un), 

has u = {ui,U2) — {u,v), X. — xi — t, n — rii — n, and Pi = pi = — 1, 
P2 = P2 = 0, p3 = p3 = 1. 

3.1 Algorithm of the tanh-method for DDEs 



Step 1 (Transform the DDE into a nonlinear DDE in T). We seek solutions 
in the traveling frame of reference, 

Q N 

= X] + X] ^j^j + S ^ d ■ n + c ■ x + S, (18) 



(17) 



where the coefficients ci, C2, . . . , cn, di, d2, ■ ■ ■ ,dQ and the phase 5 are all con- 
stants. The dot (•) denotes the Euclidean inner product. 

Using the properties of hyperbolic tangent, — tanh^n, repeatedly applying 
the chain rule, 

d^^d^dTnd^^ 2xJ^ /.qx 

dxj dxjdi^dT^ ""'^ "^UtJ ^ ' 

transforms (16) into 
A(U„+p,(T„), . . . , U„+p,(T„), U:,+p^(T„), . . . , U:,+p^(T„), 

■ ■■,u!:ip,(r„),...,u!:ip,(T„)) = o. (20) 

It is important to note that for any s {s = 1, • • • ,k), Un+p^ is a function of 
Tn and not Tn+p^ . Using the identity, 

, , , tanh X + tanh y , , 

tanh(a; + y = r^. (21 

1 + tanh a; tanh y 



we can write 



where 



_ Tn + tanh03 

1 + T„tanh0,' ^^^^ 



0s = ■ d = 'p.xdx + Vsid^ ^ V Psgdg, (23) 

and Psj is the j-th component of shift vector p^. 

Step 2 (Determine the degree of the polynomial solutions). Seeking solutions 
of the form 

U,,4T^) = Y^a,,Ti, (24) 

j=0 

we must first compute the leading exponents Mj. As in the continuous case, 
we can do this by substituting only the leading term. 



Ui^n+ps (^n) -^n+ps 



Tn + tanh 0^ 
1 + T„ tanh 0^ 



(25) 



with 0s in (23). 



Suppose we are interested in balancing terms with shift p^, then terms with 
shifts other than p^, say p^, will not effect the balance since f/j^n+p^ can be 
interpreted as being of degree zero in Tn+p,. For instance, if p; = 0, then 

^i,n+o(Tn) = T^^ is of degree M, in T„ and t/i,n+p.(T) 
degree zero in T„. 



Tn+tanh c 



1+Tn tanh <; 



Mi 

is of 
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Therefore, if we need to balance terms with shift p;, we substitute 



S = 1, 



(26) 



into (20) and proceed as in Step 2 of Section 2. We then continue with the 
union of the solutions found ior I — 1,2, . . . ,k. 

Step 3 (Derive the algebraic system for the coefficients aij). Substitute 



Mi ^ Mi 

Ui^n+ps (^n) = ^ (^ijTn+ps — ^ij 



j=0 



j=0 



Tn + tanh (\. 
1 + Tn tanh ( 



(27) 



into (20), with 0^ in (23). Applying (21) one can split tanh^s into powers of 
tanhdj. While doing so, we repeatedly clear the denominators. The resulting 
nonlinear algebraic system for the unknowns is parameterized by Ci, C2, . . . , 
cjv, tanh((ii), tanh((i2), . . . , tanh((iQ) and any parameters (lower-case Greek 
letters) in (16). 



Step 4 (Solve the nonhnear parameterized system). This step is the same as 
in the continuous case; we solve the system for in terms of the parameters 
ci, C2, . . . , Cjv, tanh((ii), tanh(o?2), ■ ■ ■ , tanh(o?Q) and any parameters (lower-case 
Greek letters) in (16). 

Step 5 (Build and test the solitary wave solutions). Substitute the solution 
found in Step 4 into (24) and reverse Step 1. Then, test the solutions (in the 
original variables) by substitution into (16). 



3.2 Example of a differential- difference equation 



To illustrate the method, we derive an exact solution of the (2-|-l)-dimensional 
Toda lattice [19], 

0^ = exp (y„_i - y„) - exp (y„ - y„+i) , (28) 

where Unix, t) is the displacement from equilibrium of the n-th unit mass under 
an exponential decaying interaction force between nearest neighbors. To write 
(28) as a polynomial DDE, set 

Oil 

= exp(y„_i - y„) - 1. (29) 
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Then, 
and (28) becomes 



exp(y„_i - y„) 



du„ 
~dt 



+ 1 



dun du 



n+1 



dt 



dt 



dt 



dxdt dt 

Integrating (31) and ignoring the integration constant, we find 

dyn 



(30) 



(31) 



dx 



(32) 



Difi^erentiating (29) with respect to x and using (30) and (32), we compute 



d 



dxdt dx 

( dU n 



(exp(y„_i - y„) - 1) 

dyn-i dyn 



dt 



dx dx J ^ 



(33) 

(34) 
(35) 
(36) 



So, (28) can be replaced by the polynomial equation (36). 
By repeatedly applying the chain rule (19) to (36), one gets 



c^c,{l-T^) 2TnK-{l-Tl)U: 



+ 



1 + ci(l - Tl)U'^ [Un-i - 2Un + = 0. (37) 
where T„ = tanh(o?in + Cix + C2t + 5). 

For this system, we have three shifts pi — —l,p2 = 0, and ps ~ 1. Substituting 
(26) into (37) and pulling off the highest exponents, we find {Mi} for shift pi 
and p3, and {Mi, Mi + 1, Mi + 2, 2Mi + 1} for shift p2. With only one term, 
neither shifts pi or p^ contribute any sohitions. Equating the two highest terms 
from the shift p2, Mi + 2 and 2Mi + 1, we find Mi = 1. 

Substituting (27) into (37), clearing the denominator and setting coefficients 
of power terms in T„ to zero, gives 

On - ci 0, 
C1C2 — tanh^(cii) — aiiC2 tanh^((ii) = 0, (38) 
C1C2 — tanh^(di) — 2aiiC2 tanh^(di) + C1C2 tanh^((ii) = 0. 
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Assuming di,ci, and an nonzero, the solution of (38) is 



aio — arbitrary, 




(39) 



On = ci = 



Then, the closed form solution of (36) is 



Un{x, t) = aio H sinh^((ii) tanh din + 




X + C2t + S 



(40) 



where aiQ,di,C2 and 5 are arbitrary. The algorithm must be repeated if any 
of the parameters in (16) are set to zero. 

3.3 Further examples 

We applied the tanh algorithm to solve a variety of nonlinear lattice equations 
in (1 + l)-dimensions. The results are summarized in Table 2. For notational 
simplicity, we denote by m„. 

The Ablowitz-Ladik equation [1,2] is a discretization of the nonlinear Schrodinger 
equation. The Toda lattices [32,34,40] describe vibrations in mass-spring lat- 
tices with an exponential interaction force. The Volterra type equations [3,17,34] 
are discretizations of the Korteweg-de Vries (KdV) and modified KdV equa- 
tions. 



4 Description of DDESpecialSolutions.m 

The format of DDESpecialSolutions is similar to the Mathematica function 
DSolve. The output is a hst of sub-lists with solutions and constraints. The 
Backus-Naur form of the function is: 

{Main Function) —> DDESpecialSolutions [(E^gwaiions), {Functions) , 

{Discrete VariaMes) , {Continuous Variables), 
{Parameters), {Options)] 
{Options) ^Verbose {Bool) |lnputForm ^ {Bool) \ 

DegreeOf ThePolynomial {List of Rules) \ 
NumericTest — > {Bool) \ SymbolicTest — > {Bool) 
{Bool) True | False 
{List of Rules) — > {m[l] — > Integer, m[2] — > Integer, ...} 
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Equation Name 


Equations 




Ablowitz-Ladik 
Lattice 

[1,2] 


Unit) = (a + UnVn)iUn+l + Un-l) - 2q;U„ 
Unit) = -(a + UnVniv„+i + Vn-l) + 2aVn 


^ «smh idi) \din + 2atsmh^idi) + 6]) 

021 

Vnit) = a2i(±l + tanh [din + 2Q;sinh^(o?i)t + 6] ) 


Toda Lattice 
[40] 


Unit) = (Un + l)('"n-l - 2Un + Un+l) 


Unit) = aio ± sinh(di) tanh[din it sinh(c/i)t + 6] 


(2+1) Dimensional 
Toda Lattice 

[19] 


i) = {^ + 1) K-l - 2Un + Un+l) 


Unix, t) = aio + ^ sinh^((ii) tanh 


'-'2 




Another Toda 
Lattice 

[32] 


Unit) = UniVn - Vn-l) 
Unit) = VniUn+1 - Un) 


Un{t) = — coth(ai jci + ci tann(ai77, + cit + d) 
Vnit) = — coth((ii)ci — ci tanh(dira + cit + S) 


Relativistic Toda 
Lattice 


Unit) = (1 + aUn)(t'„ - -Un-l) 

Vnit) = VniUn+1 - Un + aVn+1 - aVn-l) 


Unit) = ci coth(di) + ci tanh((iin + cif + 5) 

a 

Vnit) = — coth(di) tanh(din + cit + S) 

n n 


Another Relativistic 
Toda Lattice 

[32] 


Unit) = iUn+1 - Vn)Vn - - Vn-l)Vn-l 
Vnit) = VniUn+1 - Un) 


Unit) = coth(di)ci + ci tanh((iin + cit + S) 
Vnit) = coth(di)ci + ci tanh(c/in + cit + 6) 


Volterra Lattice 
[34] 


Unit) = UniVn - Vn-l) 
Vnit) = VniUn+1 - Un) 


Unit) = —ci coth((ii) + ci tanh((iin + cit + 5) 
Vnit) = —ci coth((ii) — ci tanh(din + cit + 5) 


Discretized mKdV 
Lattice 

[3] 


Unit) = (a - ul)iUn+l - Un-l) 


Unit) = =b-v/atanh(di) tanh[din + 2atanh((ii)t + S] 


Hybrid Lattice 
[17] 


Unit) = (1 + aUn + f3u'i)iUn-l - Un+l) 


u„(i)= ""^^"'2/l''^""''^'^)tanh[ 


din + tanh(di)t + 6 





Table 2 

Examples of DDEs and their solutions computed with DDESpecialSolutions.m. 
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The main function, DDESpecialSolutions, uses five subroutines corresponding 
to tlie five steps of tlie algoritlim. 

(1) DDESpecialSolutionsVarChange returns a list of nonlinear DDEs in T^. 

(2) DDESpecialSolutionsmSolver uses the result from DDESpecialSolutions- 
VarChange to determine the degree(s), Mj, of the polynomial solutions. 

(3) DDESpecialSolutionsBuildSystem uses the output of DDESpecialSolu- 
tionsVarChange and DDESpecialSolutionsmSolver to build the nonlinear 
algebraic system for the . The output of DDESpecialSolutionsBuildSys- 
tem is a list of sub-lists containing the algebraic system, the unknowns, 
the wave parameters (cij, Cj), the parameters (a, /?,... ), and the nonzero 
variables (aj,M,, di, Cj, a, P, . . .). 

(4) Algebra 'AnalyzeAndSolve AnalyzeAndSolve solves the system returned by 
DDESpecialSolutionsBuildSystem. 

(5) DDESpecialSolutionsBuildSolutions builds and tests the solutions using 
the output from Algebra AnalyzeAndSolve' AnalyzeAndSolve, DDESpe- 
cialSolutionsmSolver, and the original system of DDEs. 



5 Discussion and Conclusions 

We presented a straightforward algorithm to compute special solutions of non- 
linear polynomial DDEs, without using explicit integration. We designed the 
symbolic package DDESpecialSolutions. m to find traveling wave solutions of 
nonlinear DDEs involving the tanh functions. Our code is designed to han- 
dle DDEs with multiple continuous variables and/or multiple lattice points 
{n,m, ■ ■ ■). However, we were unable to find semi-discrete lattices with multi- 
ple lattice points in the literature. 

While the software reproduces the known solutions for many equations, there 
is no guarantee that the code will compute the complete solution set of all 
polynomial solutions involving the tanh function, especially when the DDEs 
have parameters. This is due to restrictions on the form of the solutions and 
the limitations of the algebraic solver. Furthermore, the nonhnear constraints 
which arise in solving the nonlinear algebraic system may be quintic or of 
higher degree, and therefore unsolvable in analytic form. Also, since our soft- 
ware package is fully automated, it may not return the solutions in the simplest 
form. 

Often, the nonlinear solver returns constraints on the wave parameters Cj and 
the external parameters (lower-case Greek letters), if any. In principle, one 
should verify whether or not such constraints affect the results of the previous 
steps in the algorithm. In particular, one should verify the consistency with 
the results from Step 2 of the algorithm. 
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There are several ways our algorithm could be generalized. One could look 
for solutions other than polynomials in tanh. Rational solutions would be 
an obvious choice since T^+pj is a rational expression in T^. One may try 
to compute exact solutions in terms of functions / (other than tanh). We 
restricted our algorithm to tanh solutions because /(^n+pj) = tanh(^n+pj) 
can be expressed in /(^n) = tanh(^n) via 

f{t ^- tanh(^n)+tanh(0j) _ /(^n)+tanh(0j) 

- tann^^^,; - ^^^^^^^^^^^ tanh(0,) " ! + /({„) tanh(0,) ' ^ ^ 

which allowed us to use chainrulc (19). If tanh were replaced by another func- 
tion /, the function must have the property that /(^n+p^) can be expressed as 
some function of /(Cn)- Otherwise a modified chain rule for / would no longer 
work. 
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Appendix A. Test run description 



We illustrate the use of the package DDESpecialSolutions.m on a PC. Users 
should have access to Mathematica v. 3.0 or higher. 

Put the package in a directory, say myDirectory, on drive C. Start a Mathe- 
matica notebook session and execute the commands: 

SetDirectory ["c : WmyDirectory"] ; (* specify directory *) 
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«DDESpecialSolutions .m 



(* read in the package *) 



DDESpecialSolutions [ 

{D[u[n, t] , t] == (1 + alpha*u[n, t])*(v[n, t] - v[n - 1, t] ) , 

D[v[n, t] , t] == v[n, t]*(u[n+l, t] - u[n, t] + alpha*v[n+l, t] 
alpha*v [n - 1 , t] ) 
}, {u, v}, {n}, {t}, {alpha}. 

Verbose -> False, NumericTest -> True, SymbolicTest -> True 

(*, DegreeOfThePolynomial -> {m[l] -> 1, m[2] -> 1} *)] ; 

If the DegreeOfThePolynomial {m[l] l,m[2] — > 1} were specified, the 
code would continue with this case only and not attempt to compute the 
degrees of the polynomials U{Tn) and V{Tn). 

The output of the above is: 

DDESpecialSolutionsmSolver:: "remove" : 

The potential solutions {{m[l] 1, m[2] — > 2}} are being removed because 
they are (i) negative, (ii) contain freedom, (iii) fail to balance highest expo- 
nent terms from two different terms in the original system. If Mj < 0, then 
the transformation u ^ 1/v may result in a system that DDESpecialSolu- 
tions can solve. 

Numerically testing the solutions. 

Symbolically testing the solutions. 

f 1 + a c(2) coth(c(l)) - a c(2) tanh(phase + n c(l) + t c(2)) ' 



{{{u(n,t) 



a 



c(2) (coth(c(l)) -tanh(phase + nc(l) + tc(2))) 
v[n,t) - III 

a 



The package DDESpecialSolutions. m has been tested on various PCs with 
Mathematica Yersions 3.0, 4.0, 4.1 and 5.0. 
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E(iua.ti()n Naino 


E(iua.ti()iis 


S()luti()ii(s) 


Ref. 


Ablowitz-Ladik 
Lattice 


Un{t) = {a + UnVn){UnH+Un-l)-2aUn 
Vn{t) ^-{a + UnVn{VnH+Vn-l)+2aVn 


, , asinh^fdi) / , _ 
Un{t) = ^ Mil tanh 

0-21 ^ 

Vn{t) — a2i(±l+tanh din+2a 


din+2at sinh^((ii) + 5 
smh\di)t+S]) 


) 


[1,2] 


Toda Lattice 




Un{t) — ttio ± sinh(di) tanh[din ± sinh(di)i+5] 


[40] 


(2+1) Dimensional 
Toda Lattice 


^{X,t) = {^ + 1} {Un-l-2Un + UnH) 


Un{x, t) — aio + ^ sinh^(di) tanh 


[din+''^'''^'^'^x+C2t+5 




[19] 


Another 
Toda Lattice 


Un{t) = Un{Vn-Vn-l) 
ii (fi^n 111 11 — 11 1 


Unit) = — coth((ii)ci+Ci tanh((iin+Cit + (5) 
Vn{t) — — coth((ii)ci — ci tanh(din+cii+(5) 


[32] 


Relativistic 
Toda Lattice 


11 if] — 1 1 -1- r\it 117? — 7? -1 ] 


Unit) — — — Ci coth((ii) + ci tanh((iin+Cit + 5) 
a 

Vnit) — — coth(di) tanh(din+cit + 5) 


[34] 


Another Relativistic 
ioaa Liaxtice 


Un{t) = {Urm-Vn)Vn-[Un-l-Vr^i)Vr^i 


Un(t) — coth((ii)ci + Ci tanh((iin + Cit + 5) 
i^nit) — coth((ii)ci + Ci tanh((iin+Cit+5) 


[32] 


Volterra Lattice 


U„{t) = Un{Vn-Vn-i) 
Vn{t) = Vn{Ur^i-Un) 


Unit) — —ci coth(di) +ci tanh(din+cit +5) 
Vn{t) — —ci coth(di) — ci teinh{din+cit+S) 


[34] 


Discretized mKdV 
Lattice 


Un{t) = {a-ul){Ur,^i-Un-l) 


Un{t) — ±-v/atanh(di) tanh[din+2ci;tanh(di)i+5] 


[3] 


Hybrid Lattice 




Unit) = -±V"^^^^tanh(dO ^^^^ 


Un+°'';^^ teinh{di)t+6 




[17] 



Table 3. Alternate version of Table 2. 



